RMS calculation

library(tidyverse) 
library(plotly) # for interactive plotting

Aim

This notebook demonstrates the calculation of RMS using the SAE guideline.

Calculation

SAE RMS

The SAE RMS formula definition is: \[ i_{(rms)} = \sqrt{\frac{\int_0^T i(t)^2 dt}{T}} \] Where \[i\] is the current, and \[T\] is the size of the window over which the RMS is calculated.

Pulse-Width Modulation apparent thorughput

Assuming an square signal in a canonical form, the apparent average of a waveform in a cycle T is: \[ i_{(average)} = D * i_{max} \]

Where D is the Duty Cycle calculated as the percentage of the time when the waveform is on. Since the waveform is in a canonical squared function, D is simply calculated as: \[ D = \frac{t_{on}}{t_{cycle}} \] Please note that The average of a signal is lower than the RMS. For sinusoidal signals, there is a form factor that helps to go from the average to the RMS:

Example

Signal definition

Initializing a variable with the pre-defined parameters characteristics

# Parameters
i_amp <- 18   # Instant current amplitude
pulseWidth_PW <- 50 * 1e-06 # parameter in microseconds, resulting value in seconds
cycleDuration <- 20 * 1e-03 # parameter in in milliseconds, resulting value in seconds
electrodeArea <- 1.44  # Area of the electrode in cm^2
samplingRate <- 100000  # in Hz: 100 000 data points every 1 second
pulseTrainDuration_PT <- 2 # in seconds

# Calculated parameters
delayBetweenPulses_PD <- cycleDuration - (2 * pulseWidth_PW) # in seconds
numberPulsesInTrain <- pulseTrainDuration_PT / cycleDuration # it must be an integer result for simplicity

# Wave
currentWave <- rep(
                     c(rep(i_amp, pulseWidth_PW*samplingRate),
                     rep(-i_amp, pulseWidth_PW*samplingRate),
                     rep(0, delayBetweenPulses_PD*samplingRate)),
                   numberPulsesInTrain
                  )

# Timeline
sampleNumber <- c(1:length(currentWave))

# Convert as data frame
currentWave_df <- as.data.frame(sampleNumber)
currentWave_df$currentWave <- currentWave

# Plot wave to verify shape and parameters
plot <- ggplot(data=currentWave_df, aes(x=sampleNumber, y=currentWave)) +
  geom_line(color="blue")
plot <- ggplotly(plot) # Turn to an interactive plot
plot # display the plot

Calculating the RMS per one cycle

Extracting one cycle and plotting it:

oneCycleCurrent <- currentWave_df[1:(cycleDuration * samplingRate),]

plot <- ggplot(data=oneCycleCurrent, aes(x=sampleNumber, y=currentWave)) +
  geom_line(color="red")
plot <- ggplotly(plot) # Turn to an interactive plot
plot # display the plot

Using SAE:

SAE_rms_oneCycleCurrent <- sqrt(sum(oneCycleCurrent$currentWave^2) / (cycleDuration*samplingRate))
SAE_rms_oneCycleCurrent

## [1] 1.272792

# Which is the same as:
SAE_rms_oneCycleCurrent <- sqrt(mean(oneCycleCurrent$currentWave^2))
SAE_rms_oneCycleCurrent

## [1] 1.272792

Checking manually, given our knowledge of the waveform:

# Squaring after calculating the area:
sqrt( ((2*pulseWidth_PW*i_amp)^2) / (cycleDuration) )

## [1] 0.01272792

# Squaring sample by sample
sqrt( ((2*pulseWidth_PW*i_amp^2)) / (cycleDuration) )

## [1] 1.272792

So it seems it is better not to use the area intuition, because it will lead to a different value.

Using the duty cycle definition

First, we need to put the wave into a canonical form by taking its absolute value:

oneCycleCurrent_abs <- abs(oneCycleCurrent)

# Plot to check the step
plot <- ggplot(data=oneCycleCurrent_abs, aes(x=sampleNumber, y=currentWave)) +
  geom_line(color="red")
plot <- ggplotly(plot) # Turn to an interactive plot
plot # display the plot

Calculation according to the aforementioned formula:

dutyCycle_D <- (2*pulseWidth_PW) / cycleDuration
PWM_dutyCycleAverageCurrent <- dutyCycle_D * i_amp
PWM_dutyCycleAverageCurrent

## [1] 0.09

Checking the transformation from average to RMS:

PWM_dutyCycleRMSCurrent <- (1/sqrt(dutyCycle_D)) * PWM_dutyCycleAverageCurrent
PWM_dutyCycleRMSCurrent

## [1] 1.272792

Surface area

Finally, the RMS per cycle can be simply divided by the electrode area, with either of the values

PWM_dutyCycleRMSCurrent_per_area <- PWM_dutyCycleRMSCurrent / electrodeArea
PWM_dutyCycleRMSCurrent_per_area

## [1] 0.8838835

SAE_rms_oneCycleCurrent_per_area <- SAE_rms_oneCycleCurrent / electrodeArea
SAE_rms_oneCycleCurrent_per_area

## [1] 0.8838835

This is the power per cm2. When expressing this in Watts per m^2, it is referred-to as Surface Power Density

electrodeArea_m2 <- electrodeArea / (100^2)
SAE_rms_oneCycleCurrent_per_area_m2 <- SAE_rms_oneCycleCurrent / electrodeArea_m2
SAE_rms_oneCycleCurrent_per_area_m2

## [1] 8838.835

Note regarding safety and the whole pulse train duration

The above calculations were made only per one cycle of the signal. If more than one pulse injected to the participant’s, the RMS for the whole pulse train should be calculated. If the pulse shape does not change from one cycle to the other, the RMS remains constant during the whole stimulation:

pulseTrain_TotalRMS <- sqrt(mean(currentWave_df$currentWave^2))
pulseTrain_TotalRMS

## [1] 1.272792

Comparison of the formulas to O’connor

O’Connor 2020 used the following formula to calculate the RMS:

\[ i_{rms} = I \sqrt{2 f \tau} \] Where I is the amplitude of the current in mA, and the RMS current is also in mA. f is the pulse frequency in Hz, and tau is the phase duration in seconds.

According to the examples they give on pp. 6:

# Example 1
currentAmplitude_I <- 100 # in ms
pulseFrequency_f <- 50 # in Hz
phaseDuration_tau <- 500 * 1e-6 # in s
currentRMS_I <- currentAmplitude_I * sqrt(2*pulseFrequency_f*phaseDuration_tau)
currentRMS_I

## [1] 22.36068

# Example 2
currentAmplitude_I <- 100 # in ms
pulseFrequency_f <- 20 # in Hz
phaseDuration_tau <- 500 * 1e-6 # in s
currentRMS_I <- currentAmplitude_I * sqrt(2*pulseFrequency_f*phaseDuration_tau)
currentRMS_I

## [1] 14.14214

In their paper description they say that the answers should be 22.3 mA and 14.1 mA respectively, which coincides with the calculation above.

####The section below needs to be checked according to O’Connor’s email answer

Now testing the same formula with the parameters used for the SAE and the PWM calculations:

# With the parameters above
currentAmplitude_I <- 18 # in ms
pulseFrequency_f <- 1 / (20 * 1e-03) # in Hz
phaseDuration_tau <- 20 * 1e-03 # in s
currentRMS_I <- currentAmplitude_I * sqrt(2*pulseFrequency_f*phaseDuration_tau)
currentRMS_I

## [1] 25.45584

The result seems higher, but since they do not consider an off period in between each biphasic pulse, the RMS might be overshooting.